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Pseudoscalar Meson in Two Flavors QCD with the Optimal Domain- Wall Fermion 
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We perform hybrid Monte Carlo (HMC) simulatons of two flavors QCD with the optimal domain- 
wall fermion (ODWF) on the 16"^ x 32 lattice (with lattice spacing a ~ 0.1 fm), for eight sea-quark 
masses corresponding to pion masses in the range 230-580 MeV. We calculate the mass and the 
decay constant of the pseudoscalar meson, and compare our data with the chiral perturbation 
theory (ChPT). We find that our data is in good agreement with the sea-quark mass dependence 
predicted by the next-to-leading order (NLO) ChPT, and provides a determination of the low-energy 
constants and I4, the pion decay constant, the chiral condensate, and the average up and down 
quark mass. 



Lattice QCD with exact chiral symmetry [l|, 0| is an 
ideal theoretical framework to study the nonperturbative 
physics from the first principles of QCD. However, it is 
rather nontrivial to perform Monte Carlo simulation such 
that the chiral symmetry is preserved at a high precision 
and all topological sectors are sampled ergodically. 

Since 2009, we have been using a GPU cluster (cur- 
rently constituting of 250 Nvidia CPUs) to simulate 
unquenched lattice QCD with the optimal domain-wall 
fermion (ODWF) [1, i]. Mathematically, ODWF is a 
theoretical framework which preserves the chiral symme- 
try optimally with a set of analytical weights, {ujs,s = 
1 ; • • • , iVs } , one for each layer in the fifth dimension ^ . 
Thus the artifacts due to the chiral symmetry breaking 
with finite Ng can be reduced to the minimum, especially 
in the chiral regime. The 4-dimensional effective Dirac 
operator of massless ODWF is 



D = mo[l + l^SoptiH^)], 



Sopt{Hw) — 



1+n 



T, 



1 — UJsHu 

1 + uj,,H„ 



which is exactly equal to the Zolotarev optimal rational 
approximation of the overlap Dirac operator. That is, 
Sopt{Hw) — H^RziHw), where Rz{Hw) is the optimal 
rational approximation of {H^)~^^'^ [1, @1- 

Recently we have demonstrated that it is feasible to 
perform a large-scale unquenched QCD simulation which 
not only preserves the chiral symmetry to a good preci- 
sion, but also samples all topological sectors ergodically 
To recap, we perform HMC simulations of 2 flavors 
QCD on a 16^ x 32 lattice, with ODWF at = 16 
and plaquette gauge action at /3 = 5.95. Then we com- 
pute the low-lying eigenmodes of the overlap Dirac oper- 
ator, and use its index to obtain the topological charge 
of each gauge configuration, and from which we com- 
pute the topological susceptibility for 8 sea-quark masses, 
each of 300 configurations. Our result of the topologi- 
cal susceptibility agrees with the sea-quark mass depen- 
dence predicted by the NLO ChPT and provides the 



first determination of both the pion decay constant and 
the chiral condensate simultaneously from the topologi- 
cal susceptibility. In this paper, we use the same set of 
gauge configurations to compute the mass and the decay 
constant of the pseudoscalar meson, and compare our 
results with the NLO ChPT. We find that our results 
are in good agreement with the sea-quark mass depen- 
dence predicted by NLO ChPT Q, and from which we 
obtain the low-energy constants F, E, ^3 and I4. With 
the low-energy constants, we determine the average up 
and down quark mass m^f{2 GeV), and the chiral con- 
densate I]^^(2 GeV). 

First, we outline our HMC simulation of 2 flavors QCD 
with ODWF. Starting from the ODWF action S = 
3 on the 5D lattice, we separate the even and the odd 
sites (the so-called even-odd preconditioning) on the 4D 
lattice, and rewrite V as 
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where rriq denotes the bare quark mass, denotes the 
standard Wilson Dirac operator plus a negative param- 
eter —Too (Here toq = 1.3 in this paper.), and D^^^'^^ 
denotes the part of D^^ with gauge links pointing from 
odd/even sites to even/odd sites, and 



(4 - Too) + a;-i/2(i _ + L)-^uj-^/^ 



M5 = 

L = P,L,+ P_L 



, P± = (l±75)/2, L_ 

5s-i,s', l<s<Ns 
-{mq/2mo)SN,,s' , s = 1 

Si = M5W-l/^ ^2 = (1 + L)-^uj'^/^, 
C = 1 - M.D^^MsD^^. 



Since det V = det Si^ ■ det C ■ det S2^ , and Si and 5*2 do 
not depend on the gauge field, we can just use C for the 
HMC simulation. After including the Pauli-Villars fields 
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(with mq = 2r7io)i the pseudo-fermion action for 2 flavors 
QCD (to,; ~ m^i) can be written as 

Spf = c^^Cly{CC^)-^Cpvcl). Gpv = C(2mo). (1) 

In the HMC simulation we first generate random 
noise vector ^ with Gaussian distribution, then we obtain 
(f) — CpyCS, using the conjugate gradient (CG). With 
fixed (fi, the system is evolved under a fictituous Hamilto- 
nian dynamics, the so-called molecular dynamics (MD). 
In the MD, we use the Omelyan integrator (TT|, and the 
Sexton- Weingarten multiple-time scale method 12] . The 
most time-consuming part in the MD is to compute the 
vector T] = {CC^)^^Cpv4i with CG, which is required 
for the evaluation of the fermion force in the equation of 
motion for the conjugate momentum of the gauge field. 
Here we take advantage of the remarkable floating-point 
capability of the Nvidia GPU, and perform the CG with 
mixed precision [l^. Moreover, the computations of the 
gauge force and the fermion force, and the update of the 
gauge field are also ported to the GPU. In other words, 
almost the entire HMC simulation is performed within a 
single GPU. 

Furthermore, we introduce an auxiliary heavy fermion 
field with mass mn {rriq <Si mn < 2mo), similar to the 
case of the Wilson fermion [14]. For two flavors QCD, 
the pseudofermion action (with C'h = C{mH)) becomes, 

Spf = ^^Cl,{CC^)-^CH<l> + cl>lcly{CHCl,)-^Cpvcj>H, 

which gives exactly the same fermion determinant of 
Nevertheless, the presence of the heavy fermion field 
plays a crucial role in reducing the light fermion force 
and its fluctuation, thus diminishes the change of the 
Hamiltonian in the MD trajactory, and enhances the ac- 
ceptance rate. A detailed description of our HMC simu- 
lations will be presented in a forthcoming paper jl5i] . 

Simulations are carried out for two flavors QCD 
on a 16^ X 32 lattice, for eight sea-quark masses 
m,a = 0.01,0.02,0.03,0.04,0.05,0.06,0.07, and 0.08 re- 
spectively. For the quark part, we use ODWF with 
Ng = 16, and Xmin/^max = 0.02/6.4. For the gluon part, 
we use the plaquette action at /3 = 5.95 and /3 = 5.90 
respectively. For the /3 = 5.95 ensemble, our result of 
the topological susceptibility has been presented in Ref. 

In this paper, we focus on the mass and the decay 
constant of the pseudoscalar meson for the /3 = 5.95 en- 
semble. 

For each sea-quark mass, we perform simulations on 
30 GPUs independently, with each GPU generating 400 
trajectories (Currently, the statistics has been increasing 
to 500-600 trajectories.) After discarding 300 trajecto- 
ries for thermalization, each GPU yields 100 trajectories. 
Thus, with 30 GPUs running independently, we accu- 
mulated total 3000 trajectories for each sea-quark mass. 
From the saturation of the binning error of the plaque- 
tte, as well as the evolution of the topological charge, 
we estimate the autocorrelation time to be around 10 
trajectories. Thus we sample one configuration every 10 



trajectories, then we have 300 configurations for each sea- 
quark mass. With a GPU cluster of 250 GPUs, we can 
simulate 8 sea-quark masses concurrently. It takes about 
5 months to complete the simulations for the /3 = 5.95 
ensemble. 

We determine the lattice spacing by heavy quark po- 
tential with Sommer parameter tq = 0.49 fm. The lat- 
tice spacing versus the quark mass is plotted in Fig. [T] 
Using the linear fit, we obtain the lattice spacing in 
the chiral limit, a — 0.1032(2) fm, which gives a"^ — 
1.911(4)(6) GeV, where the systematic error is estimated 
with the uncertainty of tq. 
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FIG. 1: The lattice spacing a[fm] versus mgO for two flavors 
QCD with ODWF. 

For the computation of the valence quark propagator of 
the 4D effective Dirac operator of ODWF, two methods 
have been presented in Ref. [l^. For ODWF with two 
auxiliary layers at s = and s — + 1 (with ujq — 
^N^+i — 0), and the quark fields defined in terms of these 
two auxiliary layers, then the valence quark propagator 
can be written as 

where 'D'~ly ^, is the propagator of the ODWF operator, 
which can be obtained by conjugate gradient (CG) with 
even-odd preconditioning. On the other hand, if one does 
not introduce the two auxiliary layers at s = and s = 
Ns + 1 with ujo = ujn^^i — 0, then one can obtain the 
valence quark propagator by solving the following linear 
system (with even-odd preconditioned CG), 

V{mq)\Y) = X'(2mo)B"^|source vector) (2) 

where B~]..^, ,,, = Sx,x'iP~Ss,s' + P+Ss+i,s') with periodic 
boundary conditions in the fifth dimension. Then the 
solution of ^ gives the valence quark propagator 

{Dc + mq)^^^, = (2too - mqy^ [(BY)a:,l-x',l - Sx,x'] ■ 



3 



In this paper, we compute the valence quark propaga- 
tor with the point source at the origin, and with param- 
eters exactly the same as those of the sea-quarks. 

To measure the chiral symmetry breaking due to finite 
Ng, we compute the residual mass 



.J2x (9(a;)75g(a;)g(0)75(7(0)) / 

,tr[pj +?7l,)(i?e + m,)]on, 



(3) 



{U} 



where 



denotes the chiral density at the central layer in the 5-th 
dimension, {Dc + rriq)^^ denotes the valence quark prop- 
agator with rriq equal to the sea-quark mass, tr denotes 
the trace running over the color and Dirac indices, and 
the subscript {[/} denotes averaging over an ensemble 
of gauge configurations. It turns out that, after averag- 
ing over an ensemble of a few hundreds of independent 
gauge configurations, rrires seems to be insensitive to the 
location of the origin x'^ — (0,0,0,0). Thus ([3]) gives a 
reliable estimate of the chiral symmetry breaking due to 
finite Ng. The derivation of ([3]) will be given in a forth- 
coming paper [17]. 

In Fig. [21 we plot the residual mass versus the sea 
quark mass. Using the linear fit, we obtain the residual 
mass in the chiral limit, rrires^ — 0.00040(4), less than 
5% of the lightest sea quark mass. Thus we confirm that 
the chiral symmetry is preserved to a good precision in 
our simulation. In the following, it is understood that 
each bare sea-quark mass is corrected by its residual 
mass, i.e., mo — > ma + m,res- 




FIG. 2: The residual mass versus the sea quark mass for two 
flavors QCD with ODWF. 

Using the valence quark propgator with quark mass 
equal to the sea-quark mass, we compute the time- 
correlation function of the pseudoscalar interpolator 



C^t) 



where the trace runs over the Dirac and color space. 
Then the ensemble average (C^(i)) of each mq is fitted 
to the formula 



Z 



to extract the pion mass M^ra and the decay constant 



(4) 



Fj^a — m„a 



V2Z 



where the excited states have been neglected in 

In Fig. [21 we plot M^/mg and versus m^ respec- 
tively. Here we have made the correction for the finite 
volume effect using the estimate within ChPT calculated 
up to 0{M^ / {AttFt^)*) |J^], since our simulation is done 
on a finite volume lattice with M^^L ^ 2.0 for the light- 
est sea quark, and its finite volume effect cannot be ne- 
glected. 

Taking into account of the correlation between M^/mq 
and for the same sea-quark mass, we fit our data to 
the formulas of NLO ChPT |9l 
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(5) 
(6) 



where A3 and A4 are related to the low energy constants 
Z3 and li as follows. 



h = In 



A2 



, /4 = In 



A^ 

m + 



,± = 0.140 GeV. 



The strategy of our data fitting is to search for the 
values of the parameters E, F, A3 and A4 such that they 
minimize 



1=1 



({Ml/mq), 



{Ml/mq)\ 



where Ci is the 2x2 covariance matrix for / mq and F^ 
with the same sea-quark mass, and the matrix elements 
of Ci are estimated using the binning method followed 
by the jackknife. 

For eight sea-quark masses corresponding to pion 
masses in the range 230 — 580 MeV, our fit gives 



S = [0.21855(75)(50) GeV]^ 
F = 0.08339(35)(38) GeV, 
[3 = 4. 149(35) (14), 
U = 4.582(17)(20), 



(7) 



where the systematic errors are estimated by varying the 
number of data points from 8 (M^ < 580 MeV) to 5 
(M^ < 480 MeV). 

For completeness, we re-plot everything in Fig. [31 as a 
function of M^, as shown in Fig. [H 
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FIG. 3: Physical results of 2 flavors QCD with ODWF (a) M^/m,, and (b) F^. The solid lines are the simultaneous fits to the 
NLO ChPT, for eight sea-quark masses. 
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FIG. 4: Re-plot of Fig. [3] with respect to M^. 



With the fitted parameters, we can use the NLO ChPT 
formula ([S]) to solve for the physical bare quark mass, 
with the physical inputs M^r ~ 0.135 GeV. On the other 
hand, we can also use ^ to solve for the physical bare 
quark mass, with the physical input — 0.093 GeV. 
However, the solutions of these two equations are differ- 
ent. Thus, instead of adopting either one of these solu- 
tions, we use the physical ratio 

fM^Y''^" 0.135 GcV 

— — 1.45 

J 0.093 GeV 

as the input, and solve the equation 

to obtain the physical bare quark mass m^y^ = 
0.00505(13) GeV. From ® and ©, we obtain the pion 



decay constant and the pion mass at the physical point, 

= 0.090(4) (2) GeV, (8) 
= 0.130(5)(3) GeV. (9) 

Since we have used the physical ratio 1.45 as the input, 
in principle, we can only regard either ([8|) or ([9]) as our 
predicted physical result. 

In order to convert the chiral condensate S and the 
average m„ and rud to those in the MS scheme, we cal- 
culate the renormalization factor Z^^{2 GeV) using the 
non-perturbative renormalization technique through the 
RI/MOM scheme [ll|, and our resuh is ^ 

Zf^{2 GeV) = 1.244(18)(39). 
Then the values of S and the average of rriu and nid are 
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transcribed to 

£^^(2 GeV) = [235(8)(4) MeV]^ (10) 
m^(2 GeV) = 4.06(10)(12) MeV. (11) 

Since our calculation is done at a single lattice spacing, 
the discretization error cannot be quantified reliably, but 
we do not expect much larger error because our lattice 
action is free from 0{a) discretization effects. 

We also investigated to what extent our results of the 
low-energy constants depending on the chiral symmetry 
of the valence quark propagators. We repeated above 
analysis with valence quark propagators computed with 
Ng = 32 and Xmin/^max = 0.01/6.4, which has the resid- 
ual mass mresd — 0.000191(12) in the chiral limit. The 
low-energy constants turn out to be in good agreement 
with those in ([7]). 

Moreover, our present results of the chiral condensate 
(ITUl) and the pion decay constant ^ are in good agree- 
ment with our recent results extracted from the topolog- 
ical susceptibility Q, 

=92(12)(2) MeV, 
S]^(2 GeV) = [259(6)(7) MeV]^ 

To compare our results with those obtained by other 
lattice groups, we rely on the recent review pll |. In gen- 
eral, our results of the SU{2) low-energy constants, the 
chiral condensate and the average up and down quark 
mass are compatible with those obtained by other lattice 
groups. 



Currently, for the 16^x32x16 lattice, we are increasing 
the statistics of both /? = 5.95 and f3 = 5.90 ensembles, 
targeting at a total of 8,000-10,000 trajectories for each 
sea-quark mass in each ensemble. The increased statis- 
tics will reduce the statistical error. Furthermore, we 
have been performing simulations on two larger lattices, 
20^ X 40 X 16 and 24^ x 48 x 16, for /3 = 5.95 and /3 = 5.90 
respectively. This will allow us to study the volume de- 
pendence as well as the discretization error. 

To conclude, our results of the mass and the decay con- 
stant of the pseudoscalar meson are in good agreement 
with the sea-quark mass dependence predicted by the 
next-to-leading order (NLO) ChPT, and provide a deter- 
mination of the low-energy constants I3 and Z4, the pion 
decay constant, the chiral condensate, and the average up 
and down quark mass. Together with our recent result of 
the topological susceptibility 0, these suggest that the 
nonperturbative chiral dynamics of the sea quarks are 
well under control in our HMC simulations. Moreover, 
this study also shows that it is feasible to perform large- 
scale simulations of unquenched lattice QCD, which not 
only preserve the chiral symmetry to a high precision, 
but also sample all topological sectors ergodically. This 
provides a new strategy to tackle QCD nonperturbatively 
from the first principles. 
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